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We present a framework of an auxiliary field quantum Monte Carlo (QMC) method 
for multi-orbital Hubbard models. Our formulation can be applied to a Hamiltonian 
which includes terms for on-site Coulomb interaction for both intra- and inter-orbitals, 
intra-site exchange interaction and energy differences between orbitals. Based on our 
framework, we point out possible ways to investigate various phase transitions such as 
metal-insulator, magnetic and orbital order-disorder transitions without the minus sign 
problem. As an application, a two-band model is investigated by the projection QMC 
method and the ground state properties of this model are presented. 

KEYWORDS: auxiliary field quantum Monte Carlo, multi-orbital Hubbard model, minus sign prob- 
lem, metal-insulator transition, magnetic ordering, orbital ordering 

Recently, much attention has been directed to the relevance of orbital degrees of freedom 
in strongly correlated electron systems. In addition to the coupling to the lattice distortion, 
fluctuations of the orbital component may cause fruitful interplay with those of spin degrees 
of freedom. This competition together with the strong electron correlation results in a variety 
of remarkable physical properties, for example, metal-insulator transition, spin and orbital 
order- disorder transitions and their fluctuations. 

In spite of their relevance to many aspects of physical properties in d and / electron 
compounds, theoretical as well as experimental approaches to these complex phenomena 
have not been fully developed so far. In particular, as compared with relatively simple and 
symmetric couplings in the spin degrees of freedom, the orbital degeneracy includes rather 
complicated elements which should be taken into account in theoretical models. 

One of the minimal models for this approach is the multi-orbital Hubbard model. Numer- 
ous investigations based on various approximations have been performed on this model. 
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0iSBS&&&l0S0@@EI@0@) However, for study of the serious competition men- 
tioned above and discussion of the critical properties of these phase transitions in this system, 
a more strict treatment beyond these approximations is desired. 

In this letter, we propose a framework of an auxiliary field QMC technique for this type of 
multi-orbital Hubbard models. Our technique makes it possible to study systems with any 
orbital degeneracies in any dimensions. Based on our algorithm, we identify some parameter 
regions in which no minus sign problem occurs. Using these parameter sets, we point out 
that detailed investigations on various phase transitions, such as metal-insulator transitions, 
magnetic and orbital orderings become tractable. In order to demonstrate the efficiency 
of our framework, we use it to investigate a two-band model containing level splitting and 
next-nearest neighbor hopping. Using the projection QMC method completely free from the 
minus sign problem, we discuss ground state properties of this two-band model. 

Here, we consider the multi-orbital Hubbard Hamiltonian. The explicit form of our Hamil- 
tonian is given by 

H = H t + H u + Hj + H £} (1) 

where 

w t = EEE tif ( c L<va + h.c.) (2) 

i,j v,v' a 

T~tu — X! X! (1 ~ <W<W) U vu /n iU(T ni U i a i (3) 

i U<.u' <J<.(7' 

i v<v' acr' 

+c\ yu c iula ,c\ ua ,c ivla ) (4) 

K e = Y^l e v n iv ( 5 ) 

i v 

This model has N$ sites on a bipartite lattice and Njj orbitals per site. The operator c\ va 
creates the <r-spin electron at site i — 1, • - •, N$ and in orbital v — 1, • • -, No- Here, t\- is 
the hopping integral between orbital v at site i and orbital v' at site j; U vv i is the on-site 
Coulomb interaction between orbitals v and z/; J vv i is the intra-site exchange interaction; 
and e v is the level energy for orbital v. 

Two-body interactions U uu > and J vv i are not originally independent. There is a relation 
U w = U vt £ V i—2J vv i y because of the rotational symmetry of the Coulomb terms.E10 However, 
in our model, we treat them as independent parameters. Moreover, we take U vv = U u ^ u i = 
U > and J uu > > for the convenience of our QMC algorithm. Then, we can factorize the 
two-body interaction terms (3) and (4) into the quadratic forms, as 

U ^ ( N D \ 2 U ' 
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ft J = E E -^-^iw 

i v<v' 

-(^-i)eE 

\ / i v<v' 



J vv 1 



-m. 



(7) 



Here, A% vv i — X/o-^i^'cr — So-('-'j;, -Cji/o- "I - Cj^/g-Cj^o-). 

Now, we develop the auxiliary field QMC algorithm for model (|1]). The prescription 
explained below is for the ground-state which is called the projection QMC algorithm.0'& 
It may be straightforwardly extended to the finite-temperature algorithm.il B 1 

In the ground-state algorithm, we need to calculate the density matrix, p(r; <fi) = 
(4>\ exp(— r7^)|0), where \<f>) is a trial wave function non-orthogonal to the ground state. Af- 
ter the Suzuki- Trotter decomposition of exp(— rTi) into the M = r/Ar slices, we replace the 
two-body interaction terms with non-interacting ones by introducing summations over the 
Hubbard- Stratonovich variables. The general formula of the discrete Hubbard-Stratonovich 
transformation we use here is given by 



exp [-Ar9f 2 
E E ^exp(^v^/)+0(Ar 4 



(8) 



where 7^ = 1 + ^l,rji = y2(3 — y/6l) and a = ArO > 0. Using this formula, we decouple 
the two-body interaction terms in (6) and (7), as 



n 



E E x exp 

il=±lsi=±l ^ 



iSirj h ^foii rii 



Nd 
2 



-0(Ai 



(9) 



expf-ArEE^l 



i v<v' 



nn 

i v<v' 



E E -7 iex P{^2^ 2 V / ^2^w} 



=±1 S 2 =±l 



+0(Ar 4 ), 



(10) 



where a\ = AtU/2 and ct2 = AtJ vv i/2. 

We apply these decompositions to each Suzuki- Trotter slice. Then, the one-body expres- 
sion is given as the product of these terms together with the non-interacting terms of the 
Hamiltonian (|l|), as 

p(r;0)= J2 W l ({W 2 s 2 };t;0) 

{hs 1 l 2 s 2 } 
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W l ({/isi/ 2 s 2 };r;0) 



(11) 



where W a is given by 
with 



The product of eqs. (9) and (10) over all slices amounts to the systematic error of 0(At 3 ) 
which is negligible because it is higher order than the other systematic errors from the 
Suzuki- Trotter decomposition. 

The summations over {/1S1/2S2} in eq. (11) are statistically sampled by the Monte Carlo 
technique. The product Wj-Wj, plays the role of the weight for each sample in the QMC 
updates. In general, this weight can have a negative value, which leads to the minus sign 
problem in the QMC calculation. The difficulty caused by the minus sign depends on the 
electron density, the long-range hopping and the dimensionality of the system. The details 
of this problem for this model will be reported elsewhere. 

However, the above formulation provides us with a useful property in the QMC updates. 
That is, the minus sign problem does not appear for some regions of parameter values in 
our Hamiltonian. For the parameters explained below, it is easily shown that the weight of 
the QMC sample for the up-spin is simply the complex conjugate of that for the down-spin; 
w^w l = \W^\ 2 > 0. 

The simplest case of the parameter sets which are free from the minus sign problem are 
given by the following conditions: (i) The system is at half filling; n e = N D . (ii) All 
orbitals are completely degenerated; e u = 0. (iii) The hopping integrals are non-zero only 
between the nearest neighbor sites. Under these constraints, we can easily show the complex 
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conjugate relation between the QMC weights for the up- and down-spin by the particle-hole 
transformation 

c* T - (-I)Ut (17) 
Qui -> Cj^- (18) 

This is a straightforward extension of that for the one-band Hubbard model at half filling. 

Using the symmetry about the orbital indices, the above conditions for the absence of the 
minus sign problem may easily be relaxed and generalized to the following nontrivial ones: 
(i) The system is at half filling; n e = Njj. (ii) The energy levels split symmetrically around 
the Fermi energy; e u = —e p , where v = — v + 1. (iii) The hopping integrals satisfy a 
special condition; tiff = (— l)' i ~ jf 't^ v , where \i — j\ is the Manhattan distance between site 
i and j. (iv) The intra-site exchange couplings also satisfy a special relation; J vv > = 5 Dv 'J v . 
For these sets, the particle-hole transformation to confirm the positivity of the QMC weights 
is 

Ci „ T (-iy4 T (19) 

Ciui -> c iul . (20) 

Within the former parameter sets, the perfect nesting property of the non-interacting 
energy dispersion remains unchanged, although we can change the bandwidths of each orbital 
independently. Therefore, the system is expected to become the Mott insulator for any finite 
values of U and J. Compared with this, however, the latter constraints have much more 
generality with an advantage for our purpose to investigate various phase transitions in this 
complex system. Both conditions (ii) and (iii) offer us the possibility to investigate metal- 
insulator transitions, since each of them breaks the perfect nesting property in different 
ways. These transitions are also expected to trigger order-disorder transitions of the spin 
and orbital components. 

To sum up our statement: Our QMC framework explained above provides ways of investi- 
gating various phase transitions in an iVo-orbital Hubbard model ([I]) at half filling completely 
free from the minus sign problem. Control parameters are the Coulomb interaction U, the 
intra-site exchange couplings J w i, the level splitting e v and the hopping integrals tiff with 
some constraints on the latter two to keep the particle-hole symmetry. These wide param- 
eter choices allow us to study the metal-insulator transition, the spin and orbital ordering 
transitions. 

In the following, we present an example of the applications. Here, we consider a two-band 
model with the latter choice of the above constraints; (i) n e = 2, (ii) e v= \ = —e v= 2 = £ and 
(iii) tij 1 = —Suv't for the nearest neighbor sites and tiff = —5 vv i{—l) u t' for the next-nearest 
neighbor sites. 
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This model has two trivial limits: (a) When U ^> J 3> t and t', the system is the Mott 
insulator. In this limit, the Hamiltonian is rewritten as the S — 1 spin system by the 
second-order perturbation, (b) When e ^> U,J,t and t' , since all the electrons fully occupy 
the orbital v = 2, the system becomes a band insulator. 

For finite values of U and J, the system is expected to be in the Mott insulating state at 
e = t' = because of the perfect nesting. From a consideration for the weak correlation 
limit, a finite e causes level splitting and a finite t' results in deformation of the Fermi 
surfaces of each orbital. Both should lead to self-doping between two bands. In the Hartree- 
Fock calculation, we find that these effects cause metal-insulator transitions. Depending 
on the choice of the route of the transitions, magnetic or orbital orderings take place either 
accompanied by the metal-insulator transition or independently. Therefore, the rich interplay 
of these transitions deserves study in a more reliable way. 

We present here preliminary results for the ground state properties of this model in one 
dimension (ID) calculated by the projection QMC technique to demonstrate the efficiency 
of our method. In all the calculations, we fix U/t = 2 and J/t = 0.5. We set At = t/20 
which ensures within the statistical errors the convergence to the limit At — > for all 
physical quantities we calculated. The following data are obtained from 2000 ~ 6000 QMC 
averages for the state exp(— ~H)\<f>), where we take r = 10 ~ 40t for each situation to obtain 
converged values in the ground state. The unrestricted Hartree-Fock ground state is used 
for |0) to accelerate the convergence.© 

In the following, we focus on the correlation functions for the spin and iso-spin degrees of 
freedom; 

° a (?) = ^ E<°f • Of) exp (iq ■ f l3 ) , (21) 

S hi 

where a = x, y, z. O = S and T are the spin and iso-spin operators, defined as 

^EEC'cLw (22) 

■^i = o ^ E T vu' c iv<7 C if 'oi (23) 
a vv' 

where f is the Pauli matrix. Since the Hamiltonian (P has rotational symmetry about the 
spin, we calculate the spin correlation using S = (S x + S y + S z )/3. However, since the 
iso-spin has an easy axis in the z-direction, we focus on T z . 

First, we show typical e dependence of two correlation functions at t' = in Fig. 1. For 
small values of e, the spin correlation is enhanced at q = it and the iso-spin correlation 
function has no characteristic structure compared with the non-interacting case. When e 
increases, the enhancement of S(q = 7r) disappears and the q = component of the iso-spin 
correlation grows markedly. Moreover, both correlation functions have cusp-like structures 
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at q — 2kp, where kp is the Fermi wave number of the non-interacting counterparts of these 
systems. This suggests that the level splitting leads to the self-doping from the orbital v — 1 
to 2 accompanied by the collapse of the nesting property away from e — 0. 

Next, we show the dependence on t' at e = in Fig. 2. As in Fig. 1, they also suggest a 
drastic change due to the self-doping between two orbitals. 

These QMC data clearly show the self-doping effects by the level-splitting or the next- 
nearest neighbor hopping as expected in this model. Moreover, the momentum distribution 
function clearly shows an evidence for the self-doping. We can calculate the electron density 
in each orbital by integrating the momentum distribution function for each orbital, and find 
the difference of densities between two orbitals for large values of e and t' . A more detailed 
study and analyses on the nature of these substantial changes in ID as well as in 2D will be 
reported elsewhere. 

To summarize, we have developed a new framework of an auxiliary field QMC technique for 
multi-orbital Hubbard models. Our Hamiltonian includes terms for intra- and inter-orbital 
Coulomb interaction, intra-site exchange interaction and level differences between orbitals. 
Systems with some parameter sets are completely free from the minus sign problem. For these 
parameter sets, we have pointed out the possibility to investigate various phase transitions. 
An example of the two-band model has been studied, and the ground state properties of this 
model are calculated based on our algorithm. The QMC data in one dimension clearly show 
the substantial changes induced by the self-doping between two orbitals. 

The authors thank N. Furukawa for fruitful discussions. This work is supported by a 
Grant-in-Aid for Scientific Research on the Priority Area 'Anomalous Metallic State near 
the Mott Transition' from the Ministry of Education, Science, Culture an Sports, Japan. 
The computations in this work were performed using the facilities of the Supercomputer 
Center, Institute for Solid State Physics, University of Tokyo. 
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Fig. 1. QMC results for the spin (upper-half) and iso-spin (lower- half) correlation functions at t' = for 
the two-band model in one dimension (10 sites). Each symbol denotes the difference of e in units of t as 
indicated in the inset. 
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Fig. 2. QMC results for the spin (upper-half) and iso-spin (lower-half) correlation functions at e = for 
the two-band model in one dimension (14 sites). Each symbol denotes the difference of t' in units of t as 
indicated in the inset. 
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